
import ISS_postprocessing
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt

>>>> using CPU
Running test snippet to check if MKL-DNN working
see https://pytorch.org/docs/stable/backends.html?highlight=mkl
** MKL version working - CPU version is sped up. **
>>>> using CPU
Running test snippet to check if MKL-DNN working
see https://pytorch.org/docs/stable/backends.html?highlight=mkl
** MKL version working - CPU version is sped up. **
DIPlib -- a quantitative image analysis library
Version 3.2.0 (Feb  8 2022)
For more information see https://diplib.org


samples = [
        '/media/christoffer/hfsc_processing_2/HFSC_R5/',
        '/media/christoffer/hfsc_processing_2/HFSC_R6/',
        '/media/christoffer/hfsc_processing_2/HFSC_R7/',
        '/media/christoffer/hfsc_processing_2/HFSC_R8/'
          ]


for sample in samples: 
    ISS_postprocessing.segmentation.stardist_segmentation(image_path = sample + '/preprocessing/stitched/Round4_4.tif', 
                                                      output_path = sample, 
                                                      model_name = '2D_versatile_fluo',
                                                      expand_cells = True, 
                                                      n_tiles = (4,4), 
                                                      expanded_distance = 20,)

2022-06-14 16:15:37.294989: I tensorflow/core/util/util.cc:169] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
2022-06-14 16:15:37.298974: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcudart.so.11.0'; dlerror: libcudart.so.11.0: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /home/christoffer/anaconda3/envs/3_ISS_postprocessing_new/lib/python3.8/site-packages/cv2/../../lib64:
2022-06-14 16:15:37.298992: I tensorflow/stream_executor/cuda/cudart_stub.cc:29] Ignore above cudart dlerror if you do not have a GPU set up on your machine.

Found model '2D_versatile_fluo' for 'StarDist2D'.
Loading network weights from 'weights_best.h5'.

2022-06-14 16:15:38.588043: W tensorflow/stream_executor/platform/default/dso_loader.cc:64] Could not load dynamic library 'libcuda.so.1'; dlerror: libcuda.so.1: cannot open shared object file: No such file or directory; LD_LIBRARY_PATH: /home/christoffer/anaconda3/envs/3_ISS_postprocessing_new/lib/python3.8/site-packages/cv2/../../lib64:
2022-06-14 16:15:38.588067: W tensorflow/stream_executor/cuda/cuda_driver.cc:269] failed call to cuInit: UNKNOWN ERROR (303)
2022-06-14 16:15:38.588097: I tensorflow/stream_executor/cuda/cuda_diagnostics.cc:156] kernel driver does not appear to be running on this host (borka): /proc/driver/nvidia/version does not exist
2022-06-14 16:15:38.588323: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations:  AVX2 AVX512F AVX512_VNNI FMA
To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags.

Loading thresholds from 'thresholds.json'.
Using default values: prob_thresh=0.479071, nms_thresh=0.3.
(7701, 28291)
normalize image
predict instance
1/1 [==============================] - 0s 319ms/step
1/1 [==============================] - 0s 82ms/step

  0%|                                                                                                                                                                                          | 0/16 [00:00<?, ?it/s]

1/1 [==============================] - 6s 6s/step

  6%|███████████▏                                                                                                                                                                      | 1/16 [00:06<01:32,  6.17s/it]

1/1 [==============================] - 4s 4s/step

 12%|██████████████████████▎                                                                                                                                                           | 2/16 [00:11<01:16,  5.45s/it]

1/1 [==============================] - 5s 5s/step


for sample_path in samples:  
        ad = ISS_postprocessing.annotated_objects.create_anndata_obj(spots_file = sample_path+'/decoded.csv', 
            segmentation_mask = sample_path+'stardist_segmentation_expanded.npz',#'cell_segmentation/cellpose_segmentation_expanded_2.npz' 
            output_file = sample_path+'/anndata_stardist.h5ad',
            filter_data=False, 
            metric = 'distance', 
            write_h5ad = True,
            value=  0.4,
            convert_coords = True, 
            conversion_factor = 1)

reading spots file
load coo file
assign spots to cells

AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([2, 3], dtype='int64'), …
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([2, 3], dtype='int64'), …

write h5ad
reading spots file
load coo file

AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), …
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), …

assign spots to cells
write h5ad
reading spots file
load coo file

AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([1, 2], dtype='int64'), …
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([1, 2], dtype='int64'), …

assign spots to cells
write h5ad
reading spots file
load coo file

AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), …
AnnData expects .obs.index to contain strings, but your first indices are: Int64Index([0, 1], dtype='int64'), …

assign spots to cells
write h5ad


ad =  ISS_postprocessing.annotated_objects.concat_anndata(sample_anndata_list = samples, 
                  anndata_name = 'anndata_stardist.h5ad'
                  )
ad.obs['sample_id'] = ad.obs['sample_id'].str.split('/', expand = True)[4]

/media/christoffer/hfsc_processing_2/HFSC_R5/anndata_stardist.h5ad
/media/christoffer/hfsc_processing_2/HFSC_R6/anndata_stardist.h5ad
/media/christoffer/hfsc_processing_2/HFSC_R7/anndata_stardist.h5ad
/media/christoffer/hfsc_processing_2/HFSC_R8/anndata_stardist.h5ad

Observation names are not unique. To make them unique, call `.obs_names_make_unique`.


sc.pp.calculate_qc_metrics(ad, percent_top=None, log1p=False, inplace=True)
keep = (ad.obs['total_counts'] > 5) & (ad.obs['n_genes_by_counts'] > 3)
print(sum(keep))

75168


ad = ad[keep,:]


ad.layers["raw"] = ad.X.copy()
sc.experimental.pp.normalize_pearson_residuals(ad)
normalized_exp = np.array(ad.X)
normalized_exp[normalized_exp< 0] = 0
ad.X =  normalized_exp
sc.pp.log1p(ad)

Observation names are not unique. To make them unique, call `.obs_names_make_unique`.
Observation names are not unique. To make them unique, call `.obs_names_make_unique`.


sc.tl.pca(ad, svd_solver='arpack')
plt.rcParams['figure.facecolor'] = 'white'
sc.pl.pca_variance_ratio(ad, log=True)


sc.pp.neighbors(ad, n_neighbors=20, n_pcs=40)
sc.tl.umap(ad,min_dist=0.005)


for i in [2, 3]:
    if "norm_leiden_"+str(i) in ad.obs.columns:
        plt.rcdefaults()
        with plt.rc_context({'figure.figsize': (15, 10)}):
            sc.pl.umap(ad, color = ("norm_leiden_"+str(i)),s=5,add_outline=True,legend_loc='on data',legend_fontsize=20,legend_fontoutline=2,alpha = 1)
    else: 
        print('clustering @ resolution ' + str(i))
        sc.tl.leiden(ad, resolution =i, key_added = ("norm_leiden_"+str(i)))
        plt.rcdefaults()
        with plt.rc_context({'figure.figsize': (15, 10)}):
            sc.pl.umap(ad, color = ("norm_leiden_"+str(i)),s=5,add_outline=True,legend_loc='on data',legend_fontsize=20,legend_fontoutline=2)


spatial = np.array(ad.obs[['x','y']].astype('<f8'))
ad.obsm['spatial'] = spatial
ad.obsm['xy_loc'] = spatial


cluster = 'norm_leiden_2'
for sample in ad.obs['sample_id'].unique():
    print(sample)
    ad_int = ad[ad.obs['sample_id'] == sample]    
    with plt.rc_context({'figure.figsize': (25, 17.5)}):
        sc.pl.spatial(ad_int,color=cluster ,spot_size=50)

HFSC_R5

HFSC_R6

HFSC_R7

HFSC_R8


sc.tl.rank_genes_groups(ad, cluster, method='wilcoxon', key_added = "wilcoxon")
sc.pl.rank_genes_groups(ad, n_genes=30, sharey=False, key="wilcoxon")


plt.rcdefaults()
sc.pl.rank_genes_groups_matrixplot(ad, n_genes=3, key='wilcoxon', groupby=cluster, cmap = 'viridis', )


ad.obs['project'] = 'sc'


for cluster_ in sorted(ad.obs[cluster].value_counts().head(5).index.astype(int)):
    print(cluster_)

    
    genes = list(sc.get.rank_genes_groups_df(ad,group=str(cluster_), key='wilcoxon')['names'].head(4))
    genes.insert(0, cluster)
    genes.insert(1, 'sample_id')
    plt.rcdefaults()
    sc.pl.umap(ad, color = genes, cmap = 'turbo', ncols = 3, legend_loc='on data',legend_fontsize=10,legend_fontoutline=2)
    
    ISS_postprocessing.annotated_objects.plot_specific_cluster(ad,
                    clusters_to_map = cluster,
                    cluster = str(cluster_),
                    broad_cluster = 'project',
                    key='wilcoxon',
                    size=0.5,
                    number_of_marker_genes=10,
                    sample_id_column='sample_id',
                    dim_subplots=[4, 4],)
    plt.show()

0

  
PAX3 EGFLAM RSPO3 CCL21 SPRY1 SIM1 FGF9 ID4 NKX2-8 SPON1

1

  
HOXA9 GBX2 SIM1 PTPRZ1 VEGFC NOTCH2 FGF20 COL12A1 OLIG2 HOXA11

2

  
NKX6-2 SPON1 NOTCH2 COL12A1 PTPRZ1 WNT3A FGF9 PDGFD FGF20 NKX2-8

3

  
FGF7 SPON1 WNT3A ID4 EGFLAM HOXA2 RSPO3 NKX2-8 FGF9 NOTCH2

4

  
GSX2 RSPO3 EGFLAM NOTCH2 WNT3A ID4 SPON1 SIM1 HOXA2 FGF9

import packages¶

create list of sample paths¶

run segmentation¶

create anndata¶

concat anndata objects into one file¶

filter cells¶

normalization and transformation¶